//  Date:         	08/08/2018
//  task:         	health analysis
//  project:      	World Development
//  user-written: 	renvars, interflex 

capture log close
clear all
macro drop _all
set more off
version 14.2
set matsize 800

//  #0:  setup

log using Analysis_health_Kling,replace text

//  #1: load data
use ubridge_health_long.dta

//  #2: define covariates

foreach cov in c_lpop c_age c_poverty_census c_lugbara_share c_polarization c_HHI_Religion c_literacy c_education c_employed_share c_Nonagriculture_share c_Arua_dist{
	reg treat `cov'
	}

	local controls c.c_age##treat c.c_lugbara_share##treat c.c_literacy##treat c.c_Arua_dist##treat
	su `controls'
	
	local indicators Arua_dist_miss
	tab1 `indicators'

	foreach y in monitoring_index1 effort_index1 input_index1 utilization_index1{
	center `y', replace
	drop `y'
	ren c_`y' `y'
	}

//  #3: define outcomes
local health_monitoring monitoring_index c_V_dho_visit c_V_dho_calls c_V_freq_calls c_V_BE_q67_q62 c_A_insp_qtr c_A_insp_yr
local health_effort effort_index c_V_pctpresent c_V_unauthorized_absence c_RegisterBook c_V_outreach
local health_input input_index c_V_inv_so_antimalarials c_V_inv_so_ors_analysis c_A_so_antimalarial c_A_so_ors c_A_total_so_rate c_A_Funds_Received 
local health_utilization utilization_index c_A_OutPat_Att_Total c_A_total_attendence c_A_OutPat_Ref_Total c_A_OutPat_Diag_Total c_A_TI_Total  c_A_CI_Tot c_A_ANC_IPT c_A_ANC_Iron c_A_ANC_ITN c_A_ANC_Test_Syph c_A_Mat_Admit c_A_Mat_Deliver_InUnit c_A_Mat_Moth_VitA c_A_PN_Att c_A_CH_VitA_Tot c_A_CH_Deworm_Tot

replace c_V_unauthorized_absence1=1
replace c_V_pctpresent1=1 

********************************************************************************
// #4: RUN REGRESSIONS 
********************************************************************************

foreach y in `health_monitoring'  `health_effort'  `health_input' `health_utilization' { 
	capture noisily reg `y'2 treat c.`y'1##treat healthtype
	capture est sto `y'_m1
	
	capture noisily reg `y'2 treat c.`y'1##treat `controls' `indicators' healthtype
	capture est sto `y'_m2

	capture noisily reg `y'3 treat c.`y'1##treat healthtype
	capture est sto `y'_e1
	
	capture noisily reg `y'3 treat c.`y'1##treat `controls' `indicators' healthtype
	capture est sto `y'_e2
	
* drop baseline values of the DV
	capture noisily reg `y'2 treat healthtype
	capture est sto `y'_m3
	
	capture noisily reg `y'2 treat `controls' `indicators' healthtype
	capture est sto `y'_m4

	capture noisily reg `y'3 treat healthtype
	capture est sto `y'_e3
	
	capture noisily reg `y'3 treat `controls' `indicators' healthtype
	capture est sto `y'_e4
	}
	
*************************
* // #5: Short-term tbl
*************************
forval c = 1/4  { //*************************************************************
* short term tbl
local one monitoring_index_m`c' c_V_dho_visit_m`c' c_V_dho_calls_m`c' c_V_freq_calls_m`c' c_V_BE_q67_q62_m`c' c_A_insp_qtr_m`c'
local two effort_index_m`c' c_V_pctpresent_m`c' c_V_unauthorized_absence_m`c' c_RegisterBook_m`c'  c_V_outreach_m`c'	
local three input_index_m`c' c_V_inv_so_antimalarials_m`c' c_V_inv_so_ors_analysis_m`c' c_A_Funds_Received_m`c'	
local four utilization_index_m`c' c_A_OutPat_Att_Total_m`c' c_A_total_attendence_m`c' c_A_OutPat_Ref_Total_m`c' c_A_OutPat_Diag_Total_m`c' c_A_TI_Total_m`c'  c_A_CI_Tot_m`c' c_A_ANC_IPT_m`c' c_A_ANC_Iron_m`c' c_A_ANC_ITN_m`c' c_A_ANC_Test_Syph_m`c' c_A_Mat_Admit_m`c' c_A_Mat_Deliver_InUnit_m`c' c_A_Mat_Moth_VitA_m`c' c_A_PN_Att_m`c' c_A_CH_VitA_Tot_m`c' c_A_CH_Deworm_Tot_m`c'

	tempfile short
	local n = 0
	foreach i in one two three four {

		suest ``i'' , vce(cluster cluster_id)
	
		preserve

			regsave , pval ci table(e1, parentheses(stderr) format(%8.3fc)) 
			gen grp = `n'
			keep if regexm(var , ":treat_coef") ///
				| regexm(var , ":treat_stderr") ///
				| regexm(var , ":treat_pval") ///
				| regexm(var , ":treat_ci_upper") ///
				| regexm(var , ":treat_ci_lower")
			
			replace var = subinstr(var , "treat_" , "" , .)
			gen stat = substr(var , strpos(var , ":") + 1 , .)
			replace var = substr(var , 1 , strpos(var , "_m`c'_mean") - 1)
			replace var = substr(var , 3 , .) if regexm(var , "^c_")
			
			gen order = _n
			reshape wide e1 order , i(var) j(stat) string
			renvars e1* , prefix(S)
			
			if `n' > 0 {
				append using `short'
				}
			
			save `short' , replace
		
		restore
		local ++n
		}
	
******************	
* // #6: long-term tbl
******************

local one monitoring_index_e`c' c_V_dho_visit_e`c' c_V_dho_calls_e`c' c_V_freq_calls_e`c' c_V_BE_q67_q62_e`c' c_A_insp_qtr_e`c' 
local two effort_index_e`c' c_V_pctpresent_e`c' c_V_unauthorized_absence_e`c' c_RegisterBook_e`c' c_V_outreach_e`c'	
local three input_index_e`c' c_V_inv_so_antimalarials_e`c' c_V_inv_so_ors_analysis_e`c' c_A_so_antimalarial_e`c' c_A_so_ors_e`c' c_A_total_so_rate_e`c' c_A_Funds_Received_e`c'	
local four utilization_index_e`c' c_A_OutPat_Att_Total_e`c' c_A_total_attendence_e`c' c_A_OutPat_Ref_Total_e`c' c_A_OutPat_Diag_Total_e`c' c_A_TI_Total_e`c'  c_A_CI_Tot_e`c' c_A_ANC_IPT_e`c' c_A_ANC_Iron_e`c' c_A_ANC_ITN_e`c' c_A_ANC_Test_Syph_e`c' c_A_Mat_Admit_e`c' c_A_Mat_Deliver_InUnit_e`c' c_A_Mat_Moth_VitA_e`c' c_A_PN_Att_e`c' c_A_CH_VitA_Tot_e`c' c_A_CH_Deworm_Tot_e`c'

		tempfile long
		local n = 0
		foreach i in one two three four {
		
			suest ``i'' , vce(cluster cluster_id)
			
			preserve
		
				regsave , pval ci table(e1, parentheses(stderr) format(%8.3fc)) 
		
				gen grp = `n'
				keep if regexm(var , ":treat_coef") ///
					| regexm(var , ":treat_stderr") ///
					| regexm(var , ":treat_pval") ///
					| regexm(var , ":treat_ci_upper") ///
					| regexm(var , ":treat_ci_lower")
					
				replace var = subinstr(var , "treat_" , "" , .)
				gen stat = substr(var , strpos(var , ":") + 1 , .)
				replace var = substr(var , 1 , strpos(var , "_e`c'_mean") - 1)
				replace var = substr(var , 3 , .) if regexm(var , "^c_")
	
				gen order = _n
				reshape wide e1 order , i(var) j(stat) string
				renvars e1* , prefix(L)

				if `n' > 0 {
					append using `long'
					}
			
				save `long' , replace
			
			restore
			local ++n
			}

	preserve
		use `short' , replace
		merge 1:1 var using `long'
		
		foreach k in S L {
			gen `k'_ci = "[" + `k'e1ci_lower + " , " + `k'e1ci_upper + "]" 
			replace `k'_ci = "" if missing(`k'e1ci_lower)

			destring `k'e1pval , replace
			replace `k'e1pval = `k'e1pval/2 ///
				if real(`k'e1coef) >= 0 & !missing(`k'e1pval)
			replace `k'e1pval = 1-(`k'e1pval/2) ///
				if real(`k'e1coef) < 0
			tostring `k'e1pval , replace format(%12.3f) force
			replace `k'e1pval = "" if `k'e1pval == "."

			}
		drop *lower *upper
		
		tempfile combine`c'
		save `combine`c'' , replace
	restore
	
	preserve	
		use `combine`c'' , clear
		compress

		gsort grp ordercoef -_merge
		include output_replace_edit.do
		drop order* grp _merge
		order var Se1coef Se1stderr S_ci Se1pval  Le1coef Le1stderr L_ci Le1pval 
	saveold health_analysis`c', replace
	restore
	}

// make tables	
preserve
use health_analysis1, clear
list, sep(5)
restore

preserve
use health_analysis2, clear
list, sep(5)
restore

estimates clear

*************************************************************
// #7: modertator
*************************************************************
graph drop _all
set scheme plottig, permanently

* modertator: distance
*************************************************************
* FIGURE 21 in Supplementary Information
*************************************************************

foreach y in monitoring effort input utilization {
	interflex `y'_index2 treat Arua_dist `y'_index1, title("Health `y'") xlabel("Distance to Arua") dlabel("U-Bridge") ylabel("`y'")
	graph save `y'IF, replace
}

graph combine monitoringIF.gph effortIF.gph inputIF.gph utilizationIF.gph, col(2)


********************************************************************************
* Robustness I: multi-level models (school random effects)
********************************************************************************

// prepare data in long format
set more off
preserve
local health_monitoring monitoring_index* c_V_dho_visit* c_V_dho_calls* c_V_freq_calls* c_V_BE_q67_q62* c_A_insp_qtr* c_A_insp_yr*
local health_effort effort_index* c_V_pctpresent* c_V_unauthorized_absence* c_RegisterBook* c_V_outreach*
local health_input input_index* c_V_inv_so_antimalarials* c_V_inv_so_ors_analysis* c_A_so_antimalarial* c_A_so_ors* c_A_total_so_rate* c_A_Funds_Received* 
local health_utilization utilization_index* c_A_OutPat_Att_Total* c_A_total_attendence* c_A_OutPat_Ref_Total* c_A_OutPat_Diag_Total* c_A_TI_Total*  c_A_CI_Tot* c_A_ANC_IPT* c_A_ANC_Iron* c_A_ANC_ITN* c_A_ANC_Test_Syph* c_A_Mat_Admit* c_A_Mat_Deliver_InUnit* c_A_Mat_Moth_VitA* c_A_PN_Att* c_A_CH_VitA_Tot* c_A_CH_Deworm_Tot*
local covs healthtype c_age c_lugbara_share c_Arua_dist

keep `covs' `health_monitoring'  `health_effort' `health_input' `health_utilization' treat cluster_id

local health_monitoring monitoring_index c_V_dho_visit c_V_dho_calls c_V_freq_calls c_V_BE_q67_q62 c_A_insp_qtr c_A_insp_yr
local health_effort effort_index c_V_pctpresent c_V_unauthorized_absence c_RegisterBook c_V_outreach
local health_input input_index c_V_inv_so_antimalarials c_V_inv_so_ors_analysis c_A_so_antimalarial c_A_so_ors c_A_total_so_rate c_A_Funds_Received 
local health_utilization utilization_index c_A_OutPat_Att_Total c_A_total_attendence c_A_OutPat_Ref_Total c_A_OutPat_Diag_Total c_A_TI_Total  c_A_CI_Tot c_A_ANC_IPT c_A_ANC_Iron c_A_ANC_ITN c_A_ANC_Test_Syph c_A_Mat_Admit c_A_Mat_Deliver_InUnit c_A_Mat_Moth_VitA c_A_PN_Att c_A_CH_VitA_Tot c_A_CH_Deworm_Tot

reshape long `health_monitoring'  `health_effort' `health_input' `health_utilization' , i(cluster_id) j(wave)

saveold "health_long.dta", replace

// estimate

tempfile monitoring_index1 monitoring_index2 monitoring_index3 monitoring_index4 
tempfile effort_index1 effort_index2 effort_index3 effort_index4
tempfile input_index1 input_index2 input_index3 input_index4
tempfile utilization_index1 utilization_index2 utilization_index3 utilization_index4

foreach y in monitoring_index effort_index input_index utilization_index {
	
	mixed `y' i.treat healthtype || cluster_id:, r
	capture est sto `y'_m5
	
	mixed `y' i.treat healthtype i. wave || cluster_id:, r
	capture est sto `y'_m6
	
	mixed `y' i.treat healthtype if wave!=3 || cluster_id:, r
	capture est sto `y'_m7
	estadd local control "no"
	margins , dydx(treat) post
	parmest, label bmat(r(b)) vmat(r(V)) saving(``y'1')
	
	mixed `y' i.treat healthtype c.c_age##treat c.c_lugbara_share##treat c.c_Arua_dist##treat if wave!=3 || cluster_id:, r
	capture est sto `y'_m8
	estadd local control "yes"
	margins , dydx(treat) post
	parmest, label bmat(r(b)) vmat(r(V)) saving(``y'2')
	
	mixed `y' i.treat healthtype if wave!=2 || cluster_id:, r
	capture est sto `y'_m9
	estadd local control "no"
	margins , dydx(treat) post
	parmest, label bmat(r(b)) vmat(r(V)) saving(``y'3')
	
	mixed `y' i.treat healthtype c.c_age##treat c.c_lugbara_share##treat c.c_Arua_dist##treat if wave!=2 || cluster_id:, r
	capture est sto `y'_m10
	estadd local control "yes"
	margins , dydx(treat) post
	parmest, label bmat(r(b)) vmat(r(V)) saving(``y'4')
	
	}

list, sep(5) 
	
restore

preserve
use `monitoring_index1', clear 
	append using `monitoring_index2' 
	append using `monitoring_index3' 
	append using `monitoring_index4' 
	append using `effort_index1' 
	append using `effort_index2' 
	append using `effort_index3' 
	append using `effort_index4' 
	append using `input_index1' 
	append using `input_index2'
	append using `input_index3' 
	append using `input_index4'
	append using `utilization_index1'
	append using `utilization_index2'
	append using `utilization_index3'
	append using `utilization_index4'
	keep if parm=="1.treat"
	
	gen n=_n
	gen CovAdjust="No"
	replace CovAdjust="Yes" if n==2 | n==4 | n==6 | n==8 | n==10 | n==12 | n==14 | n==16
	
	gen var="Monitoring index"
	replace var="Effort index" if n==5 | n==6 | n==7 | n==8  
	replace var="Input index" if n==9 | n==10 | n==11 | n==12
	replace var="Utilization index" if n==13 | n==14 | n==15 | n==16
	
	gen year=1
	replace year=2 if n==3 | n==4 | n==7 | n==8 | n==11 | n==12 | n==15 | n==16
	drop parm label
	
saveold "health_multilevel", replace	
restore

estimates clear

*************************************************************
* * Robustness II: permutation tests
*************************************************************
set more off
use ubridge_health_long.dta, clear

foreach i in monitoring_index effort_index input_index utilization_index {
	regress `i'2 treat `i'1 healthtype
	mat beta=e(b)
	svmat double beta, names(matcol)
	ren betatreat beta`i'2
	drop beta`i'1 betahealthtype beta_cons
	su beta`i'2
	local beta`i'2 = r(mean)
	ritest treat _b[ treat ], reps(5000) right seed(05022018) saving(`i'2, replace): regress `i'2 treat `i'1 healthtype, r

	regress `i'3 treat `i'1 healthtype
	mat beta=e(b)
	svmat double beta, names(matcol)
	ren betatreat beta`i'3
	drop beta`i'1 betahealthtype beta_cons
	su beta`i'3
	local beta`i'3 = r(mean)
	ritest treat _b[ treat ], reps(5000) right seed(05022018) saving(`i'3, replace): regress `i'3 treat `i'1 healthtype, r
	
	preserve
			use `i'2, clear
			gen t= `beta`i'2'
			gen c=0
			replace c=1 if _pm_1 >= t
			egen p = mean(c)
			saveold `i'2, replace
	restore
	
	preserve
			use `i'3, clear
			gen t= `beta`i'3'
			gen c=0
			replace c=1 if _pm_1 >= t
			egen p = mean(c)
			saveold `i'3, replace
	restore
	}

****************************************		
set scheme plottig, permanently
preserve 
graph drop _all

* Short-term
foreach i in monitoring effort input utilization{
	use `i'_index2.dta, clear
	quietly su t, de
	local t=round(r(mean),.001)
	quietly su p, de
	local p=round(r(mean),.001)

		hist _pm_1 , kdensity lcolor() title(Health `i' (year-1)) addplot(pci 0 `t' 5 `t', lc(red) lp(solid) lw(thick)) text(1.5 .4 "pval=`p'", size(small)) text(4 `t' "{&tau}=`t'",place(nw) size(small)) name(`i'2) xlabel(-.5(.1).5) xsca(r(-.5 .5)) legend(off)
	}
restore 	

* long-term
preserve 

****************************************
* FIGURE 3 in Supplementary Information
****************************************	
foreach i in monitoring effort input utilization{
	use `i'_index3.dta, clear
	quietly su t, de
	local t=round(r(mean),.001)
	quietly su p, de
	local p=round(r(mean),.001)

		hist _pm_1 , kdensity lcolor() title(Health `i' (year-2)) addplot(pci 0 `t' 5 `t', lc(red) lp(solid) lw(thick)) text(1.5 .4 "pval=`p'", size(small)) text(4 `t' "{&tau}=`t'",place(nw) size(small)) name(`i'3) xlabel(-.5(.1).5) xsca(r(-.5 .5)) legend(off)
	}
restore 
	
	graph combine monitoring2 monitoring3 effort2 effort3 input2 input3 utilization2 utilization3, title("Health (unweighted indices)") col(2)
	
log close
clear
exit